A land-based controlled-source electromagnetic method for oil field exploration: An example from the Schoonebeek oil field

Controlled-source electromagnetic (CSEM) data are sensitive to the subsurface resistivity distribution, but 3D inversion results are ambiguous, and in-depth interpretation is challenging. Resolution and sensitivity analysis as well as the influence of noise on resolution have been used to quantify 3D inversion performance. Based on these numerical studies, a land-based CSEM survey was designed and carried out at the Schoonebeek oil field, the Netherlands. The acquired data were processed and subsequently inverted for the resistivity distribution. The 1D and 3D inversion of horizontal electric-field data show the reservoir at the right depth, matching well-log data without using a priori knowledge about the actual reservoir depth. We used a 1D model with fine layering as a starting model for 3D inversion. Synthetic data inversions and sensitivity tests demonstrate that resistive or conductive bodies inside the reservoir zone may be well-detectable with our limited acquisition geometry. Spatial variations in the reservoir resistivity are visible in the measured data and after inversion by assuming good knowledge of the background resistivity distribution. The reservoir resistivity and size, however, have to be interpreted with care considering the intrinsically low resolution of electromagnetic (EM) which is further reduced by manmade EM noise. INTRODUCTION Since the first survey of an oil field offshore Angola using marine controlled-source electromagnetic (CSEM) sounding for determining the absence or presence of hydrocarbons in a known reservoir (Ellingsrud et al., 2002), the CSEM method has gained a lot of interest in the industry. The marine CSEM method is generally recognized to give additional information and is especially applied for hydrocarbon exploration, primarily with the objective to derisk drilling activities (Darnet et al., 2007; Constable, 2010; Fanavoll et al., 2014). Although in marine CSEM surveys, dense 2D profiles or complex 3D grids with tens of transmitter tow lines and hundreds to thousands of receiver deployments are feasible (Constable, 2010), such source and receiver coverage is very difficult to achieve for terrestrial surveys because of various logistical constraints (access restrictions, urbanization, presence of infrastructure, and noise sources) and typically limited equipment availability. In addition, deploying sources of sufficient strengths is difficult. In contrast to the marine case in which the transmitter is situated in the most conductive area, the seawater, land transmitters are usually situated on the earth’s surface, and the source electrodes are deployed in a medium much more resistive than sea water, which limits achievable source current amplitudes. Furthermore, more near-surface heterogeneities exist on land than at the seafloor and near-electrode heterogeneities can result in more complex electric-field patterns. Although electromagnetic (EM) methods were initially developed on land (Streich, 2016), the mentioned problems and challenges of land EM make it difficult to replicate the marine EM exploration success on land. To date, the seismic method is the principal geophysical method that is routinely applied for hydrocarbon exploration on land. Seismic wave-propagation method provides higher resolution than the diffusive CSEMmethod. However, CSEMmethods have the advantage to be sensitive to resistive objects and to spatial variations in the resistivity. They provide generally higher resolution than potential-field methods such as gravity (Li and Oldenburg, 1998; Dell’Aversana et al., 2012). Large-scale EM surveys are not commonly carried out, but many wells exist with resistivity logs. Gathering additional EM data would provide knowlManuscript received by the Editor 10 January 2017; revised manuscript received 23 August 2017; published ahead of production 13 November 2017; published online 28 December 2017. Delft University of Technology, Department of Geoscience and Engineering, Delft, The Netherlands. E-mail: [email protected]; G.G.Drijkoningen@ tudelft.nl. Shell Global Solutions International BV, Amsterdam, The Netherlands. E-mail: [email protected]. GFZ German Research Centre for Geosciences, Potsdam, Germany. E-mail: [email protected]. © 2018 Society of Exploration Geophysicists. All rights reserved. WB1 GEOPHYSICS, VOL. 83, NO. 2 (MARCH-APRIL 2018); P. WB1–WB17, 26 FIGS. 10.1190/GEO2017-0022.1 D ow nl oa de d 01 /2 5/ 18 to 1 31 .1 80 .1 30 .2 42 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / edge of the resistivity distribution inside an oil or gas field, which is essential to evaluate its processes during hydrocarbon production or steam injection. Numerical feasibility studies on land exist (Wirianto et al., 2010) and show that subsurface responses are weaker compared with marine CSEM responses and are more complicated to interpret due to near-surface inhomogeneities. Up to now, only a few academic EM surveys and a small number of industry surveys were conducted, which can be considered technology trials rather than commercial surveys (Streich and Becken, 2011; Grayver et al., 2014; Tietze et al., 2015; Streich, 2016). Grayver et al. (2014) showed that for a land CSEM survey, it is possible to get good 3D inversion results consistent with regional geology. They deployed 39 5C receivers along an 11 km long line centered at a CO2 injection site and eight CSEM transmitters well distributed around the area. Several 3D inversion algorithms have been developed over the last decade to reconstruct resistivity models from marine and land EM data (Commer and Newman, 2008; Plessix and Mulder, 2008; Grayver et al., 2013; Oldenburg et al., 2013; Schwarzbach and Haber, 2013; Zhdanov et al., 2014). A recent review of the latest developments can be found in Newman (2014). In this work, we study the applicability of land CSEM for recovering the resistivity distribution inside a hydrocarbon reservoir by 1D and 3D inversions. The underlying mathematical theory of EM inversion is well-established, but the practical application to real data requires thorough analysis of uncertainties in the inversion result. Therefore, we numerically investigate resolution capability in terms of the survey geometry, noise distortion, reservoir depth, and resistivity distribution inside the reservoir. We use a deterministic inversion approach, in which gradients derived from the forwardmodel process are used to update a sequence of models. Receivers are either inverted separately or together, and their inversion results are compared. Based on the results of the analysis, a target-oriented land CSEM acquisition field setup was designed for the Schoonebeek region in the Netherlands where steam is injected into an oilbearing reservoir at a depth of approximately 700 m for enhanced oil recovery. The obtained data were inverted for the resistivity distribution inside the reservoir where the injected steam may lead to more complex resistivity patterns close to the injector and producer wells. We show how a sparse source-receiver configuration has the potential of resolving complex resistivity patterns inside the reservoir zone. Furthermore, we discuss to what extent small-scale variations of resistivity can be detected. INVERSION ALGORITHMS 1D inversion We try to find a model for which we can compute data that will best fit the measured data. Finding this best-fitting model requires iterative (forward) modeling in which after every iteration, the model can be updated. This was done using a reflectivity forward modeling code (Streich and Becken, 2011; Streich et al., 2011). To minimize the number of iterations necessary to find the best-fitting model, we use a Gauss-Newton-type method to compute model updates based on an objective function that we seek to minimize. The objective function is given by f1D 1⁄4 kWdðd − FðmÞÞk2 þ μkWmð∂zðm −m0ÞÞk; (1) where the first term is the misfit of the model’s computed response FðmÞ to the data d that are the real and imaginary parts of the measured electric-field data. The data vector d can contain data from multiple transmitters, receivers, and frequencies. The data weights are defined as Wd 1⁄4 diagðwiÞ, where the weights wi are calculated by multiplying error estimates with the absolute value of the data. They weigh the relative contribution of each datum to the misfit. Data with large errors are down weighted to limit their influence, whereas data with small errors will have a larger impact on the total misfit. The second term is a norm of the model roughness and is computed by applying a differencing operator ∂z, a matrix of first derivatives with respect to depth, to the elements of the model vector m, andm0 is a reference model, e.g., the starting model. We use a bounded logarithmic transform of the conductivity σh;v to the model parameters (Grayver et al., 2013), where σh;v represents horizontal and vertical conductivities. For the logarithmic transform, upper and lower conductivity bounds of 10−4 S∕m (10;000 Ωm) and 1.5 S∕m (2∕3 Ωm) were used. The variable Wm is a diagonal weighting matrix that weighs the model variations and can consist of different measures of the model norm (Farquharson and Oldenburg, 1998); μ is a regularization parameter that weighs the dataand model-dependent terms of the objective function and contains a reduction exponent that gradually decreases the regularization (“cooling”) at each iteration (Haber et al., 2000). The regularization parameter at iteration n is computed as μ 1⁄4 ðμiÞ∕ðnþ 1Þp where the regularization parameter μi weighs the dataand model-dependent terms and p is a parameter that defines by how much the influence of regularization is reduced at each iteration. For the modeling study, a simple least-squares weighting of the model variations μi 1⁄4 0.05 and p 1⁄4 1.67 was used. By minimizing the first-order derivative of the conductivity depth profile in addition to minimize the data misfit, the regularization seeks to generate a smooth model. Although there is no physical argument in using a smoothing constraint, smooth models are less likely to result in over-interpretation of the data because they will not contain smallscale features that are poorly constrained by the data. Our algorithm is similar to the Occam inversion (Constable et al., 1987; deGroot Hedlin and Constable, 1990; Key, 2009), but it does not include a search for the optimum regularization parameter. We found that in most cases, field data cannot be fit to within a prescribed error level and searching for an adequate regularization parameter would return very small values of the objective function, leading to instable models. The algorithm is stable and rapidly convergent: A smoothed version of the true structure is typically recovered in 12–16 iterations. Themisfit xrms is defined as the global root-mean-square (rms) error: xrms 1⁄4 1 n ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Xn i1⁄41 Wdi 1⁄2di − FiðmÞ 2 s ; (2) where n is the number of data points. The inversion is terminated if the target rms is reached or if either the objective function or the rms cannot be reduced during several subsequent iterations. 3D inversion The forward modeling algorithm used in the 3D inversion code consists of solving the second-order partial differential equations (Mulder, 2006; Plessix et al., 2007): WB2 Schaller et al. D ow nl oa de d 01 /2 5/ 18 to 1 31 .1 80 .1 30 .2 42 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / iωμ0σ̄E − ∇ × ∇ × E 1⁄4 −iωμ0Js (3) where Eðω; xÞ represents the electric-field components as a function of angular frequency ω and position x, the complex conductivity σ̄ 1⁄4 σ þ iωε includes conductivity σðxÞ and electric permittivity ε, and μ0 is the magnetic permeability; μ and ε are assumed to be constant at their free-space values. We do not consider frequency-dependent conductivity variations; induced polarization effects are irrelevant for the fairly resistive geologic settings studied here. The current source is Jsðω; xÞ. Equation 3 is discretized using a finite-volume-type discretization scheme resulting in a linear system of the form AE 1⁄4 F, where A is the discretized Maxwell operator, E represents the vector of the electric-field components on the discretized model, and F is the source vector. A conjugategradient iterative scheme, the BiCGStab2 scheme, preconditioned by a multigrid solver (Mulder, 2006) is used to solve the large and sparse, but symmetric and non-Hermitian, matrix A. For multifrequency, multisource CSEM modeling, the code is parallelized over sources and frequencies. For further details of the forward-modeling engine, the reader is referred to Plessix et al., (2007). The inverse problem seeks to find a conductivity model that minimizes the weighted least-squares functional: f3D 1⁄4 1 2 kWdðd − FðmÞÞk2 þ R: (4) The value of FðmÞ is obtained by solving equation 3, and d contains the electric-field data. To balance the update of the shallow and deeper parts, a model weighting scheme that mainly depends on depth is used (Plessix and Mulder, 2008). We invert for the logarithm of the conductivity and impose upper and lower hard conductivity bounds. The data weights Wd are computed from the data amplitudes and noise; R is a regularization term given by R 1⁄4 P αnðjm −m0jÞ, where n are the spatial directions and αn are positive numbers that are calibrated such that the magnitude of the regularization term remains a small fraction of the total objective function value throughout the inversion. The objective function is then solved by using a quasi-Newton method, the limitedmemory BFGS method (Byrd et al., 1995), and a box average filter is applied to the gradient of the objective function to smooth the spatially rapid variations and geometry imprints arising from the sparse source and receiver spacing or fine model discretization. For computational efficiency, we terminated all inversions after a maximum of 25 iterations, when in most cases, convergence had slowed down, and we were not expecting further significant model updates and misfit reductions. In the following, vertical transverse isotropic (VTI) inversions were carried out for 1D and 3D inversion using horizontal electric-field components Ex and Ey under the assumption that horizontal and vertical resistivities are different. If not mentioned explicitly, we focus on the vertical and do not show the horizontal resistivity because the inversion is more sensitive to the vertical resistivity inside the resistive reservoirs of interest. SYNTHETIC 3D RESOLUTION TESTS Resolution and sensitivity analysis The Schoonebeek reservoir layer is located at approximately 700–800 m depth, and its structure is well-known from well logs, 3D seismic data, and production data. We found from forward modeling tests that the strongest reservoir responses should be obtained, while maintaining acceptable signal amplitudes, at source-receiver offsets between approximately 3 and 5 km. The detailed numerical resolution and sensitivity analysis were carried out based on known properties of the Schoonebeek oil field in the Netherlands with the objective of obtaining high resolution toward a resistive reservoir target zone. Resolution is the measure for the distance that two separated objects should have, to be identified as two objects and not as one single object. Thus, two objects can be resolved when their separation is larger than the minimum separation distance and their contrast to the background is higher than a predefined value. An object is detectable when amplitude differences of data with and without this object are sufficiently large and the phase differences are above the detection limit. An object can therefore be detectable but not resolvable. We modeled a range of complex structures inside the Schoonebeek reservoir target zone to understand (1) the minimum size of a feature that is still detectable taking into account a specific acquisition geometry limited by various logistical constraints and (2) resolution limits and the capability to recover resistivity values. To minimize the impact and cost of the survey, we aim for a sparse field setup with few source locations and a small number of receivers. Synthetic sources and receivers were placed at positions where actual field deployment is possible such that the survey covers one of the locations where steam is being injected into the reservoir (Figure 1). The source available for our survey was a transmitter that feeds currents with a fixed phase relationship into three grounded electrodes. The overall source polarization can be adjusted by applying a constant phase shift to the three source currents (see Appendix A). We refer to this phase shift as the source polarization angle. Notice that the actual spatial orientation of the source polarization is determined by the combination of source geometry and the polarization angle applied. For the following synthetic calculations, if not mentioned otherwise, we modeled data for two transmitter positions and 15 receivers using source polarization angles of −30°, 30°, and 270°. For forward modeling tests and as a background model for our 3D studies, we use a 1D background resistivity model that was obtained by taking the shallow subsurface resistivity information for the top 150 m from regional well-log data and resistivity at greater depths from 1D real-data CSEM inversion results of the Schoonebeek region (Figure 2). We tried a range of other starting models, but they gave poor convergence of the 3D inversion for field data, and thus, we limit our study to the mentioned background model. In Figure 3, we display examples of data d1 calculated for the model shown in Figure 2, and data d2 for a model containing a